How the cheetah got his spots

Why aren't they striped?

Fig: from ChatGPT

Patterns

How do patterns form?

How patterns form (pattern formation) and how they are recognised by the brain (pattern recognition) are important areas of study in Complex Systems.

We have already seen fractals, which show self-similar structure across scales (e.g. in plant branching or river networks).

But there are many other processes in the natural world that form patterns.

Rayleigh–Bénard convection: where heat flow generates stable, symmetric roll patterns in fluids.

The pattern that forms in a hot coffee on a cold morning might look like a bit like this garlic but they are very different mechanisms. The semi-regular pattern of garlic cloves is the result of a growth process producing a such as honeycomb-like structures (Voronoi tesselation) in biological tissues or crystals.

Fig from: @alizauf (on X)

Elastic instabilities: are thought to underlie folding patterns in the cerebral cortex of higher animals

Travelling waves, fronts, oscillations, and clustering: as seen in chemical reactions, neural activity, or population dynamics

from IPython.display import YouTubeVideo
display(YouTubeVideo("jRQAndvF4sM", width=400))
# https://www.youtube.com/watch?v=jRQAndvF4sM

Diffusion-based mechanisms: including reaction–diffusion systems that generate spots, stripes, and spirals.

This is what we will focus on and in particular, one of the most celebrated mechanisms to explain self-organising spatial structures... Turing reaction-diffusion instability as a source of symmetry breaking.

Alan Mathison Turing

The chemical basis of morphogenesis

Turing was ahead of his time in many ways.

Turing, Alan M. The Chemical Basis of Morphogenesis. 1952. Philosophical Transactions of the Royal Society; 237: 37-72

Reaction-diffusion systems

Turing proposed that complex spatial structures may arise from the interaction of two types of chemical substances, 'morphogens' diffusing through tissue.

It is the interaction of two morphogens: an "activator" and an "inhibitor" with competing interests that allows for the autonomous generation of spatial patterns during development via a type of symmetry breaking.

This 'dappled' pattern in was created by Turing 'in a few hours by a manual computation':

Reaction

We use $A$ and $B$ to denote the activator and inhibitor respectively and will explore the specific example of the following chemical reactions:

$$ A+2B\rightarrow 3B $$

i.e. a single species of $A$ will interact with two species of $B$ to form three species of $B$ (sort of like predator-prey system)

This is part of the Gray-Scott model.

For this interaction to happen (with some probability) the particles first have to come into contact with each other, which happens via diffusion.

Diffusion

But first... Brownian motion (aka Wiener processes) was named for Robert Brown (1827) who noticed that pollen grains suspended in water moved in a move in a jittery, random fashion.

The underlying cause remained a mystery until Albert Einstein (1905) explained it using the kinetic theory of gases, providing a crucial bridge between the microscopic world of atoms and the macroscopic world we observe.

Once Loop Reflect

The equation Einstein derived to describe Brownian motion is:

$$ \langle x^2 \rangle = 2Dt, $$

Returning to the example we are considering of the two morphogens...

Let the diffusion of $B$ (a single particle depicted with white-blue-purple), $D_B$ to be twice that of $A$ (a single particle depicted with yellow-orange-brown), $D_A$.

Turing instabilities

Turing examined a chemical system in stable equilibrium, where, in the absence of spatial effects (variations in concentration of particles across space), no net reactions were occurring. He demonstrated that this stability could be disturbed solely by introducing diffusion.

Regulation

We will also include:

• Feed rate, $f$: controls the supply of $A$
• Kill rate, $k$: controls the decay of $B$ (converted to $C$ via the chemical reaction: $B \rightarrow C$, where $C$ is an inert product of these irreversible reactions)

This is also lumped in with the 'reaction' part of the dynamics (because it isn't diffusion).

Putting it all together

So what happens in this system?

Complexity Explorables

Returning to why a cheetah has spots

Many animals display distinct and repeatable patterns—stripes, spots, or swirls—yet they share a remarkably similar set of genes. So how does this visual diversity arise?

The key lies not in each species having a unique genetic "blueprint" for their patterns, but in how biological shape and scale influence the behaviour of how the morphogens spread and interact to stimulate or suppress development.

The beauty of mathematical modelling is that once you're confident your model captures the biology, you can start to play around with it and ask biological questions that are difficult to answer through experiments alone.

— Kit Yates

I've sort of lied to you...

Remember...

Fig from: The first paragraph of Turing's paper on morphogenesis.

Is the world discrete or continuous?

If it is discrete, when is pretending it is continuous acceptable/helpful?

If it is continuous, when is pretending it is discrete acceptable/helpful?

The true nature of reality might encompass both aspects, and our understanding is still evolving as scientific theories develop.

Is Brownian motion discrete or continuous?

Aspect	Nature	Explanation
Particle	Discrete	You’re tracking the motion of an individual microscopic particle (e.g. a pollen grain), which is a discrete object.
Space	Continuous	The particle’s position $x(t)$ is represented as a real-valued variable in continuous space.
Time	Continuous (in theory)	Although early explanations used time steps, the modern model treats time as continuous.
Motion	Continuous but non-smooth	The motion is continuous in time but highly irregular — it's modelled as a Wiener process, which is continuous everywhere but differentiable nowhere.

Modelling Brownian motion in a discrete world

Brownian motion — the jittery movement of a microscopic particle suspended in fluid — can be understood as the result of many small, random steps.

This is by a random walk.

Modelling diffusion in a continuous world

The alternative is to aggregate particles.

If we do this, then the spatio-temporal dynamics can be modelled with partial differential equations in a continuous field.

e.g., we can consider a flux of particles that is proportional to the concentration gradient and describe how that concentration changes with time as particles diffuse from regions of high concentration to regions of low concentration.

Gif from Wiki: Molecular diffusion from a microscopic and macroscopic point of view

Check your foundational understanding:
What exactly is the gradient? What exactly is the Laplacian?

```{toggle}  my answers...

The gradient ($\nabla f$) is a vector that points in the direction of steepest increase of a scalar field and gives the rate of change in that direction.
The Laplacian ($\nabla^2 f$) is a scalar that measures the net rate of flow in or out of a point — i.e. the sum of the second derivatives. It tells you if the point is a peak, valley, or flat relative to its surroundings.

In short:
- **Gradient**: “Which way is up?”
- **Laplacian**: “Am I higher or lower than my neighbours?”
```

Transport equation

$$ \frac{\partial u}{\partial t} + \nabla \cdot (\mathbf{v} u) = D \nabla^2 u + R(u) $$

• $u(\mathbf{x}, t)$: the scalar quantity being transported (e.g. concentration)
• $\mathbf{v}$: velocity field
• $D$: diffusion coefficient
• $R(u)$: source, sink or interaction term
• $\nabla \cdot (\mathbf{v} u)$: advection (transport by flow)
• $D \nabla^2 u$: diffusion (spreading out due to random motion)

Return to the reaction-diffusion system

$$ \frac{\partial U}{\partial t}=R_U(U,V) + D_U \nabla^2 U $$

$$ \frac{\partial V}{\partial t}=R_V(U,V) + D_V \nabla^2 V $$

"Solving" the PDEs

The chemical reactions and diffusion-driven pattern formation is happening at the particle level (discrete).

But... Spatio-temporal dynamics is typically described via particle concentrations with PDEs (continuous).

But... We need to actually simulate this (discrete again!).

Diffusion in a coarse-grained (discretised) world...

This doesn’t require full calculus: partial differential equations can be translated directly into code using simple update rules. e.g., the concentration fields $U$ and $V$ evolve over time as:

$$ U \rightarrow U + D_U \, \Delta U $$

$$ V \rightarrow V + D_V \, \Delta V $$

Here, $\Delta U$ and $\Delta V$ represent the change in concentration due to diffusion, often approximated numerically using the discrete Laplacian, denoted $\nabla^2 U$ and $\nabla^2 V$.

A note on choosing colour schemes

Numerical approximation of the gradient using the finite-difference method

Thus, the gradient vector at $(i,j)$ is approximated by:

$$ \nabla U(i,j) \approx \left( \frac{U(i+1,j) - U(i-1,j)}{2},\; \frac{U(i,j+1) - U(i,j-1)}{2} \right) $$

This gives the rate of change of $U$ in the $x$- and $y$-directions respectively, based on values from both sides of the point.

Numerical approximation of the Laplacian with the finite-difference method

To approximate the Laplacian $\nabla^2 U$ of a scalar field $U$ on a grid, we use the finite-difference method to compute second-order partial derivatives in both spatial directions.

$$ \nabla^2 U(i,j) \approx U(i+1,j) + U(i-1,j) + U(i,j+1) + U(i,j-1) - 4U(i,j) $$

This operation can be expressed as a convolution with the following kernel (or stencil matrix):

$$ L=\left[ \begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{array} \right] $$

This discrete Laplacian matrix serves the same purpose as the continuous Laplacian operator — quantifying local spatial variation — and is commonly used in modelling diffusion.

Regulation

In addition to diffusion, chemical reactions regulate the local concentrations of particles by adding, removing, or transforming them. These processes happen can be incorporated directly into the model as additional update rules at each time step.

feed:

$A$ particles are fed into the system. Let's assume that particles are added to each cell, at each time step.

$$ U(i,j)\rightarrow U(i,j) + f(1-U(i,j)) $$

We also remove particles $B$ (and $C$) at a rate proportional to their current concentration of the particular grid cell to keep the total density of particles fixed:

$$ V(i,j)\rightarrow V(i,j) - fV(i,j) $$

kill:

$B$ particles are converted to inert $C$ particles at a rate, $k$, proportional to their current concentration at the particular grid cell.

$$ V(i,j) \rightarrow V(i,j) - kV(i,j) $$

Reaction

The likelihood of a particular reaction occurring depends on how frequently particles collide within a given grid cell — which in turn depends on their local concentrations.

Consider:

$$ A + 2B\rightarrow 3B $$

$$ U \rightarrow U - UV^2 $$

$$ V \rightarrow V + UV^2 $$

(I’ve omitted the $(i,j)$ subscripts here for clarity, but note that these updates are applied locally at each grid cell.)

Putting it all together:

The reaction + diffusion of particles was modelled with PDEs:

$$ \frac{\partial U}{\partial t}= D_U \nabla^2 U + f(1-U)-UV^2 $$

$$ \frac{\partial V}{\partial t}=D_V \nabla^2 V - (k+f)V +UV^2 $$

that were discretised for numerical simulation:

$$ U \rightarrow U + D_U\Delta U + f(1-U) - UV^2 $$

$$ V \rightarrow V + D_V\Delta V - (k+f)V+ UV^2 $$

Simulating

Models like this can't be solved analytically.

But that's ok because we don't really want to solve it.

i.e. we don't really care where the stripes on a zebra are, only that for a particular set of parameters and conditions stripes occur.

There are some beautiful simulators out there to play with. See for example:

• Simulator 1
• Simulator 2

But we would still like to create our own... You'll do this week's Workshop.

So is the world discrete or continuous?

Modelling Brownian motion with a random walk

This is a microscopic and stochastic description of the motion of an individual particle as a series of discrete, random steps in space and time.

Aspect	Nature	Explanation
Time	Discrete	Einstein modelled Brownian motion as a sequence of random jumps in small time intervals.
Steps	Discrete	In a random walk, the particle takes discrete steps in space at each discrete time interval.
Transition	Probabilistic	Each step is random — often modelled as drawn from a normal distribution with zero mean, leading to no net drift over time.
Limit	Continuous stochastic process	As the time step and the step length scales appropriately, the random walk converges to Brownian motion, which is the basis of the diffusion equation.

Modelling Brownian motion with the diffusion equation

This is a macroscopic and deterministic description of how a distribution of particles evolves smoothly over time in continuous space.

Aspect	Nature	Explanation
Time	Continuous	Time is treated as a continuous variable, allowing smooth evolution of the system over any time scale.
Space	Continuous	Space is modelled as a continuous domain; the concentration of particles is defined at every point in space.
State	Density field	Instead of tracking individual particles, the model describes the concentration of many particles as a scalar field $u(x,t)$.
Transition	Deterministic	The diffusion equation (a partial differential equation) describes how the concentration evolves over time based on physical laws, such as Fick’s law.
Limit	Emerges from stochastic behaviour	The diffusion equation arises as the continuum limit of a large number of independent random walks, capturing the average behaviour of many particles undergoing Brownian motion.